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We apply a recently introduced method for global optimization to deter- 
mine the ground state energy and configuration for model metallic clusters. 
The global minimum for a given iV-atom cluster is found by following the 
damped dynamics of the N particle system on an evolving potential energy 
surface. In this application, the time dependent interatomic potential inter- 
polates adiabatically between the Lennard-Jones (LJ) and the Sutton-Chen 
(SC) forms. Starting with an ensemble of initial conditions corresponding 
to the ground state configuration of the Lennard-Jones cluster, the system 
asymptotically reaches the ground state of the Sutton-Chen cluster. We de- 
scribe the method and present results for specific cluster size N = 15, when 
the ground state symmetry of LJ^v and SCtv differ. 
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I. INTRODUCTION 



Determination of the lowest energy configuration for a cluster of N atoms is a nontrivial 
task ]I],^|,[5y3]^] . The complexity arises in part from the exponentially (in N) large number 
of minima in the potential energy surface ||. Furthermore, the geometry of the potential 
energy landscape itself can make it computationally hard. 

The problem is, however, simple to describe: Find the lowest energy minimum of a N 
body potential energy surface, V(r), 

V{v) = Y,V{nj) (1) 

hi 

where r are the atomic coordinates, and V(rij) is the interatomic potential of interaction 
between atoms % and j. For small N one can hope to enumerate all possible minima and 
decide the lowest of these, but even for moderate N and for the simplest V such as the 
Lennard- Jones (LJ) potential, this becomes difficult. A case in point is the 38 atom LJ 
cluster which has the so-called "double funnel" structure; the global minimum, which has 
octahedral symmetry, is marginally lower than the first excited state which has icosahedral 
symmetry. These were respectively found by the basin hopping technique J7| and a genetic 
algorithm method ||||. 

A number of techniques of global optimization have been applied to this problem |T|, p!0|JTl 



and by now there are extensive compilations of global minima for a number of different 
clusters |I2| , notably those described by two-body or many-body potentials which are 
commonly applied in atomistic simulations. A major difficulty is in ensuring that the algo- 
rithms reach the global minimum without being trapped in local minima. One method of 
overcoming such trapping |13| is by transforming the PES, broadening the thermodynamic 



transitions so as to increase the probability of finding the global minimum at temperatures 



where the free energy barriers are almost unsurmountable. For example, addition of a linear 
term to the PES provides a compressing effect which has been shown to be successful in 



locating the true minima in multiple funneled global structures |L4]. Locatelli and Schoen 



| 15fl used such transformations to locate the global minimum for non-icosahedral clusters. 
We have recently proposed a new method of global optimization wherein time-dependence 



is introduced in the potential energy landscape |IB| . The evolving landscape is designed in 
a manner such that asymptotically, the potential energy surface develops into the surface 
of interest. A number of other techniques can be used to follow the evolving minima. 
Applications have been made to determining the ground state configurations of simple cluster 
systems [|TI| . 

In the present paper we apply this method to determine the ground state configurations 
and energies of atomic clusters described by the many-body Sutton-Chen potential |T7|JT8 



by switching from a known ground state. Initially, the interaction is chosen to be the 
Lennard- Jones potential, 

Vo(?) = £ Vb(r y ) = £4e[(a/r^) 12 - (a/r^} (2) 

ij ij 

while the surface of interest is the potential 

v/(?) = eE[^E(^r-^j, "=E(£r ( 3 ) 



One choice for the time dependent potential energy surface is [JTSj 



V(t) = V (v)g(t) + V f (v)h(t) (4) 

with g(t) an adiabatically varying switching function that interpolates between 1 and 0, and 
h{t) doing the reverse. 
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In the next section we describe the method as applied to the problem of ground state 
energy determination for Sutton-Chen clusters. Detailed results are presented for one cluster 
size, while the more general application and results are indicated in brief. This is followed 
by a discussion and summary. 

It is a pleasure to dedicate this article to Steve Berry who has directly and indirectly 
influenced much of the development in the area of cluster studies over the past few decades. 
We have learned a lot from him, both in conversation as well as through his many articles 



and reviews 19 



II. ADIABATIC OPTIMIZATION 



The adiabatic optimization method [jig] , is a heuristic technique for locating minima. 
The essential idea is as follows. 

Time dependence is introduced into the potential energy landscape directly by the in- 
corporation of slowly varying terms as discussed in Eq. (Q). A given choice is made for the 
switching functions g(t) and h(t), though in practice the choice does not affect the results 
greatly. A similar application of the adiabatic principle to determine semiclassical ground 



states of multidimensional systems [^0] has noted the insensitivity of the technique to the 
precise form of the switching function, so long as the induced variation of the potential 
energy surface is slow enough. We note, parenthetically, that the switching principle has 
wide applicability, and in recent work has been used in the computation of the free-energy 
of finite clusters [plf . 

Location of the evolving minima can be done by any of a number of techniques. The 
simplest procedure is to introduce damping into the equations of motion and allow the 
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system to evolve to a position of rest in a potential minimum; by starting with an ensemble 
of initial configurations and varying the available parameters, a number of minima can 
be located, and the putative global minimum can be recognized. Elsewhere [16| we have 



suggested the conjugate gradient [^] or simulated annealing (SA) |23| as other possible 
methods for locating the minima. It is likely that of these, the conjugate gradient technique 
will be more efficient as compared to SA though some SA variants |23] may also provide a 
suitable method for following the evolving minima. 
The overall procedure can be summarized as follows: 



1. Take the initial configuration of the N atom cluster to be the ground state for the LJ 



/V 



cluster 12 . 



2. Choose some switching function, say g{t). Similarly choose h(t), and the simplest 
choice, which we make here, is h(t) = I — g(t). We have explored a large variety of 
switching functions and in the present application we use g{t) = cos 2 (37r££) exp(— £t), 
where ( is the adiabaticity parameter. 

3. Perform molecular dynamics simulations for this iV-particle cluster with forces deriving 
from Eq. (|J), with an additional damping term, namely the equations of motion 

dV 

mr k + 7 r fc + — = 0, k = 1, 2 . . . , N, (5) 

where is the position vector for the kth particle, m is its mass and 7 is the damping 
constant. 

4. Vary ( and 7, keeping in mind the natural timescales of the problem. Evolve to a 
minimum energy configuration, namely when the particle velocities become zero; the 



5 



lowest energy found in an ensemble of simulations is the ground state energy predicted 
by the present method. 



A. Results for Sutton Chen global minima 

Here we attempt to switch from the minimum of the LJ^v system to the minimum of the 
SCtv system. Both sets of minima have been extensively studied earlier and are tabulated 
in the Cambridge Cluster Database |12| . A point of interest is that for the Sutton-Chen 9-6 



family of potentials, the symmetries of the global minimum configurations are fre- 

quently different from the symmetries of the Lennard- Jones minima, so that in the adiabatic 
switching process, the cluster atoms must also move so as to adopt a different symmetry. 

We present detailed results for the cluster size N = 15, though we have applied this 
technique to larger clusters and obtained results in agreement with the current standards 
|12|1 . For the 15 atom LJ cluster, the ground state has the point group symmetry C2 V while 



for the SC cluster the symmetry is Dq^- Shown in Fig. 1 is a plot of the potential energy 
versus time for a particular choice of g{t), ( and 7. Also shown is the effect of instantaneous 
switching, namely taking the limit g(t) = 0, where it can be seen that the system finds the 
nearest available local minimum from which it does not move. The time-dependence in the 
potential effectively permits the system to explore the multidimensional potential energy 
landscape of the SC cluster in an efficient manner. Finding a local minimum does not trap 
the system since there is always kinetic energy until the adiabatic switch is essentially over. 
Inset in the figure is a schematic of the cluster configuration at different times during the 
process, showing how the cluster both contracts as well as rearranges to eventually reach 
the minimum of the SC surface. 
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FIG. 1. Plot of the potential energy versus iteration number, switching from the Lennard Jones 
potential to the Sutton Chen. The switching function used is g(t) = exp (-Ct)cos 2 (3vrCt) with 
£ = 0.4 and 7 = 0.1. The time step is 0.01 in units natural to the LJ cluster, for which we also 
take e = a = 1. At time t = the cluster has the Ci v symmetry. At different times, as indicated, 
the cluster configuration is shown, and asymptotically, the configuration reached is the 9-6 Sutton 
Chen global minimum, with Dqj symmetry. The parameters used for this latter model are taken 
from [12]. The dashed line shows the result of the simulation in the absence of switching, namely 
when the LJ potential is suddenly transformed to the Sutton Chen potential. 
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It should be added that we have performed simulations for a variety of cluster sizes and 
in all cases we find that the procedure successfully finds the tabulated minima of SC clusters; 
these are not presented here since the details are repetitive. As we have emphasized else- 



where [0, the present method is heuristic, and thus some exploration of different switching 
functions, variation in the adiabaticity and damping parameters, and indeed the choice of 
initial potential, Vq(?) is necessary. 



III. SUMMARY AND DISCUSSION 

In this paper we have presented the outline of a general procedure for global optimization 
with specific application to the problem of cluster ground state geometry determination. 
The application here, to the determination of the minimum of model metallic (Sutton- 
Chen) clusters by adiabatically deforming the potential energy surface relevant to model 
rare-gas (Lennard Jones) clusters is meant to be illustrative rather than exhaustive: the 
method introduced here is one of a class of techniques that employs time-dependence in 
the potential energy surface to enhance the exploration of phase space in contrast to other 
means of achieving the same objective ||10|| . 

A multiplicity of techniques is needed to approach hard problems such as global optimiza- 
tion. Few rigorous results are available, and application of most techniques is not guaranteed, 
with few possible exceptions, to give reliable (or certifiable) results. The present adiabatic 
switching method locally solves the optimization for an evolving surface, and thus mimics 
other methods of making large scale excursions in configuration or phase space. 

We are presently studying this technique in detail with respect to the variation of pa- 
rameters as well as to functional variations. One of the main issues of concern, and one that 
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we are addressing in current work, is the relative efficiency of this method in comparison to 
other global optimization techniques. In a number of applications, we find that this method 
gives very encouraging results, and permits the determination of fairly reliable minimum 
energy configurations for a wide variety of cluster systems. The flexibility of choice of a 
number of starting potentials including the free particle case |1| as well as the flexible 



choice of switching functions and parameters, and finally the flexibility in the dynamical 
evolution all combine to suggest that while the method is heuristic, it holds promise. 
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